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Abstract 

We address the problem of the stability of the computations of resultants and subresultants 
of polynomials defined over complete discrete valuation rings ( e.g. or k[[t]] where k is a 

field). We prove that Euclide-like algorithms are highly unstable on average and we explain, 
in many cases, how one can stabilize them without sacrifying the complexity. On the way, 
we completely determine the distribution of the valuation of the subresultants of two random 
monic p-adic polynomials having the same degree. 


1 Introduction 

As wonderfully illustrated by the success of Kedlaya-type counting points algorithms [9], p-adic 
technics are gaining nowadays more and more popularity in computer science, and more specifically 
in symbolic computation. A crucial issue when dealing with p-adics is those of stability. Indeed, 
just like real numbers, p-adic numbers are by nature infinite and thus need to be truncated in 
order to fit in the memory of a computer. The level of truncation is called the precision. Usual 
softwares implementing p-adics {e.g. MAGMA [5], pari [3], SAGE [11]) generally tracks the precision 
as follows: an individual precision is attached to any p-adic variable and this precision is updated 
after each basic arithmetic operation. This way of tracking precision can be seen as the analogue 
of the arithmetic intervals in the real setting. We refer to §2.1.2 for more details. 

In the paper [6], the authors propose a new framework to control p-adic precision. The aim 
of this paper is to illustrate the technics of loc. cit. on the concrete example of computation of 
GCDs and subresultants of p-adic polynomials. There is actually a real need to do this due to 
the combination of two reasons: on the one hand, computating GCDs is a very basic operation for 
which it cannot be acceptable to have important instability whereas, on the other hand, easy ex¬ 
perimentations show that all standard algorithms for this task {e.g. extended Euclide’s algorithm) 
are very unstable. Figure 1 illustrates the instability of the classical extended Euclide’s algorithm 
(c/ Algorithm 1) when it is called on random inputs which are monic 2-adic polynomials of fixed 
degree (see also Example 2.3). Looking at the last line, we see that extended Euclide’s algorithm 
outputs the Bezout coefficients of two monic 2-adic polynomials of degree 100 with an average loss 
of 160 significant digits by coefficient whereas a stable algorithm should only loose 3.2 digits on 
average. This “theoretical” loss is computed as the double of the valuation of the resultant. Indeed 
Cramer-like formulae imply that Bezout coefficients can be computed by performing a unique di¬ 
vision by the resultant, inducing then only the aforementioned loss of precision (see §2.1.2, Eq. (8) 


Degree 

Loss of precision (in number of significant digits) 

Euclide algorithm 

expected 

5 

6.3 

3.1 

10 

14.3 

3.2 

25 

38.9 

3.2 

50 

79.9 

3.2 

100 

160.0 

3.2 


Figure 1: Average loss of precision when computing the GCD of two random monic polynomial of 
fixed degree over Z 2 . 
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Algorithm 1 : Extended Euclide’s algorithm 

Input : Two polynomials A,Bg Q p [X] (whose coefficients are known at given precision) 
Output: A triple D , U, V such that D = AU + BV = gcd(A, B) 


1 

2 

3 

4 

5 

6 

7 

8 


Si «- A; «- 1; Vi «- 0 
S'2 i — B ; U 2 t— 0; V 2 <— 1 
fc 2 

while Sk ^ 0 do 

Q, Sk+i <— quotient and remainder in the Euclidean division of Sk- 1 by Sk 
Uk+i t— Uk~ 1 — QUk 
Vk+i — QI4 

k <r- k + 1 


9 return S k -i,U k -i,V k - 1 


21 

IE 

7T 

I< 

k 

2l<n[A'] 

2l<n[X] 

2ln[X] 

Res dA ’ dB (A,.B) 

Res^’ dB (A,S) 


a commutative ring (without any further assumption) 
a complete discrete valuation ring 
a uniformizer of IE 
the fraction field of IE 
the residue field of IE 

the free 21-module consisting of polynomials over 21 of degree < n 
the free 21-module consisting of polynomials over 21 of degree < n 
the affine space consisting of monic polynomials over 21 of degree n. 

The resultant of A and B “computed in degree (gU, c?b)” 

The j-th subresultant of A and B “computed in degree (gU,g(b)” 

Figure 2: Notations used in the paper 


for a full justification). Examining the table a bit more, we observe that the “practical” loss of 
precision due to Euclide’s algorithm seems to grow linearly with respect to the degree of the input 
polynomials whereas the “theoretical” loss seems to be independant of it. In other words, the 
instability of Euclide’s algorithm is becoming more and more critical when the degree of the input 
increases. 


Content of the paper The aim of this article is twofold. We first provide in §3 a theoretical 
study of the instability phenomenon described above and give strong evidences that the loss of 
precision grows linearly with respect to the degree of the input polynomials, as we observed empir¬ 
ically. In doing so, we determine the distribution of the valuation of the subresultants of random 
monic polynomials over Z p (c/ Theorem 3.3). This is an independant result which has its own 
interest. 

Our second goal, which is carried out in §4, is to rub out these unexpected losses of precision. 
Making slight changes to the standard subresultant pseudo-remainder sequence algorithm and 
using in an essential way the results of [6], we manage to design a stable algorithm for computing 
all subresultants of two monic polynomials over Z p (satisfying an additional assumption). This 
basically allows to stably compute GCDs assuming that the degree of the GCD is known in advance. 

Notations Figure 2 summerizes the main notations used in this paper. The definitions of many 
of them will be recalled in §2. 


2 The setting 

The aim of this section is to introduce the setting we shall work in throughout this paper (which 
is a bit more general than those considered in the introduction). 
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2.1 Complete discrete valuation rings 

Definition 2.1. A discrete valuation ring (DVR for short) is a domain W equipped with a map 
val: W — > Z U {+ 00 } — the so-called valuation — satisfying the four axioms: 

1. val(x) = +00 iff x = 0 

2. val(xy) = val(x) + val(y) 

3. val(x + y) > min(val(x), val(y)) 

4. any element of valuation 0 is invertible. 

Throughout this paper, we fix a discrete valuation ring W and assume that the valuation on 
it is normalized so that it takes the value 1. We recall that W admits a unique maximal ideal m, 
consisting of elements of positive valuation. This ideal is principal and generated by any element 
of valuation 1. Such an element is called a uniformizer. Let us fix one of them and denote it by w. 
The residue field of W is the quotient W/m = W/ttW and we shall denote it by A:. 

The valuation defines a distance d on W by letting d(x,y) = e - val O-y) f or a ]i x,y £ W. 
We say that W is complete if it is complete with respect to d, in the sense that every Cauchy 
sequence converges. Assuming that W is complete, any element x £ W can be written uniquely as 
a convergent series: 

X = X 0 + Xl7T + X 2 7T 2 H-+ x n Tr n + ■■■ (1) 

where the xf s lie in a fixed set S of representatives of classes modulo 7 r. Therefore, as an additive 
group, W is isomorphic to the set of sequences N —> k. On the contrary, the multiplicative structure 
may vary. 

Let K denote the fraction field of W. The valuation v extends uniquely to I\ by letting 
val(^) = val(x) — val(y). Moreover, it follows from axiom 4 that K is obtained from W by 
inverting n. Thus, any element of K can be uniquely written as an infinite sum: 

OO 

X = XiTT 1 (2) 

i=io 

where i 0 is some relative integer and the xfs are as above. The valuation of x can be easily read 
off this writing: it is the smallest integer i such that Xi ^ 0 (mod 7r). 

2.1.1 Examples 

A first class of examples of discrete valuation rings are rings of formal power series over a field. 
They are equipped with the standard valuation defined as follows: val(]C ?;>0 ait 1 ) is the smallest 
integer i with cq ^ 0. The corresponding distance on &[[£]] is complete. Indeed, denoting by f[i] 
the term in t l in a series f £ fe[[t]], we observe that a sequence (/ n )„>0 is Cauchy if and only if the 
sequences (f n [i])n >0 are all ultimately constant. A Cauchy sequence (f n )n >0 therefore converges 
to £V >0 a it 1 w h ere a i is the limit of /„[*] when n goes to + 00 . The DVR fc[[t]] has a distinguished 
uniformizer, namely t. Its maximal ideal is then the principal ideal (t) and its residue field is 
canonically isomorphic to k. If one chooses 7 x = t and constant polynomials as representatives 
of classes modulo t, the expansion (1) is nothing but the standard writing of a formal series. 
The fraction field of fc[[i]] is the ring of Laurent series over k and, once again, the expansion (2) 
corresponds to the usual writing of Laurent series. 

The above example is quite important because it models all complete discrete valuation rings 
of equal characteristic, i.e. whose fraction field and residue field have the same characteristic. On 
the contrary, in the mixed characteristic case (i.e. when the fraction field has characteristic 0 and 
the residue field has positive characteristic), the picture is not that simple. Nevertheless, one can 
construct several examples and, among them, the most important is certainly the ring of p-adic 
integers Z p (where p is a fixed prime number). It is defined as the projective limit of the finite 
rings Z/p”Z for n varying in N. In concrete terms, an element of Z p is a sequence (x n ) n >o with 
x n £ r Llp n 'L and x„+i = x n (mod p n ). The addition (resp. multiplication) on Z p is the usual 
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coordinate-wise addition (resp. multiplication) on the sequences. The p-adic valuation of ( x n ) n >o 
as above is defined as the smallest integer i such that Xi ^ 0. We can easily check that Z p equipped 
with the p-adic valuation satisfies the four above axioms and hence is a DVR. A uniformizer of 
Z p is p and its residue field is Z/pZ. A canonical set of representatives of classes modulo p is 

{0,1,...,p— 1}. 

Given a p-adic integer x = ( x n ) n >o , the i-th digit of x n in p-basis is well defined as soon as 
i < n and the compatibility condition £ n +i = x n (mod p n ) implies that it does not depend on 
n. As a consequence, a p-adic integer can alternatively be represented as a “number” written in 
p-basis having an infinite number of digits, that is a formal sum of the shape: 

a 0 + a\p + a 2 p 2 H-+ a n p n H- with oq G {0,1,... ,p — 1}. (3) 

Additions and multiplications can be performed on the above writing according to the rules we 
all studied at school (and therefore taking care of carries). Similarly to the equal characteristic 
case, we prove that Z p is complete with respect to the distance associated to the p-adic valuation. 
The writing (3) corresponds to the expansion (1) provided that we have chosen n = p and S = 
{0,1,...,p — 1}. The fraction field of Z p is denoted by Q p . 

2.1.2 Symbolic computations over DVR 

We now go back to a general complete discrete valuation ring W, whose fraction field is still 
denoted by K. The memory of a computer being necessarily finite, it is not possible to represent 
exhaustively all elements of W. Very often, mimicing what we do for real numbers, we choose 
to truncate the expansion (1) at some finite level. Concretely, this means that we work with 
approximations of elements of W of the form 


JV-l 

x = XiTT 1 + 0(n N ) with N G N (4) 

i—0 

where the notation 0(ir N ) means that the Xi’s with i > N are not specified. 

Remark 2.2. From a theoretical point of view, the expression (4) does not represent a single 
element x of W but an open ball in W, namely the ball of radius e~ N centered at Xiir 1 (or 

actually any element congruent to it modulo n N ). In other words, on a computer, we cannot work 
with actual p-adic numbers and we replace them by balls which are more tractable (at least, they 
can be encoded by a finite amount of information). 

The integer N appearing in Eq. (4) is the so-called absolute precision of x. The relative precision 
of x is defined as the difference N—v where v denotes the valuation of x. Continuing the comparison 
with real numbers, the relative precision corresponds to the number of significant digits since x 
can be alternatively written: 

N—v—1 

x = p v + 0(n N ) with yj = Xj+ V and z/o 7 ^ 0 . 

3 =0 


Of course, it may happen that all the Xi s (0 < i < N ) vanish, in which case the valuation of x is 
undetermined. In this particular case, the relative precision of x is undefined. 

There exist simple formulas to following precision after each single elementary computation. 
For instance, basic arithmetic operations can be handled using: 


(a + 0(n Na )) +(b + 0(n Nb )) = a + b + 0(n miniN *’ Nb) ), 

(a + 0( tt n “)) ~(b + 0(7r Nb )) = a — b + 0( 7 r min ( JV “’ JVi ’)), 

(a + 0( 7 r JV “)) x (b + 0{n Nb )) = a 6 + 0(7r min ( JV “ +val ^’ JVi > +val ( a »). 

(a + 0(7T JVa )) ^ (6 + 0{TT Nb )) = ^+ 0( 7r 'nin(^-val(6),JV 6 +val(a)-2val(b))^l 


(5) 

( 6 ) 

(7) 

( 8 ) 
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with the convention that val(a) = N a (resp. val( 6 ) = Nf,) if all known digits of a (resp. b) are zero. 
Combining these formulas, one can track the precision while executing any given algorithm. This 
is the analogue of the standard interval arithmetic over the reals. Many usual softwares (as SAGE, 
magma) implement p-adic numbers and formal series this way. We shall see later that this often 
results in overestimating the losses of precision. 

Example 2.3. As an illustration, let us examine the behaviour of the precision on the sequence 
(Ri) while executing Algorithm 1 with the input: 

A = A 5 + (27 + 0( 2 5 ))X 4 + (11 + 0(2 5 ))A 3 + (5 + 0(2 5 ))A 2 + (18 + 0( 2 5 ))X + (25 + 0(2 5 )) 
B = A 5 + (24 + 0(2 5 ))A 4 + (25 + 0(2 5 ))A 3 + (12 + 0(2 5 )) A 2 + (3 + 0(2 5 )) A + (10 + 0(2 5 )). 

The remainder in the Euclidean division of A by B is £3 = A — B. According to Eq. ( 6 ), we do 
not loose precision while performing this substraction and the result we get is: 

S 3 = (3 + 0(2 5 )) A 4 + (18 + 0(2 5 )) A 3 + (25 + 0(2 5 )) A 2 + (15 + 0(2®)) A + (15 + 0(2 5 )). 

In order to compute S 4 , we have now to perform the Euclidean division of S 2 = B by S 3 . Noting 
that the leading coefficient of S 2 has valuation 0 and using Eq (5)—( 8 ), we deduce that this operation 
does not loose precision again. We get: 

S 4 = (26 + 0(2 5 )) A 3 + (17 + 0(2 5 )) A 2 + (4 + 0(2 5 )) A + (16 + 0(2 5 )). 

We observe now that the leading coefficient of S 4 has valuation 1 . According to Eq. ( 8 ), divising 
by this coefficient — and therefore a fortioti computing the euclidean division of S 3 by S 4 — will 
result in loosing at least one digit in relative precision. The result we find is: 

S 5 = (| + 0(2 2 )) A 2 + (6 + 0(2 3 )) A + (3 + 0(2 3 )) . 

S. v .✓ S, v ✓ Si v .✓ 

rel. prec.=4 rel. prec.=2 rel. prec.=3 

Continuing this process, we obtain: 

S 6 = (20 + O(2 5 ))A+(l2 + 0(2 5 )) and S 7 = 1 + 0(2). 

The relative precision on the final result S 7 is then 3, which is less than the initial precision which 
was 5. 

2.2 Subresultants 

A first issue when dealing with numerical computations of GCDs of polynomials over W is that the 
GCD function is not continuous: it takes the value 1 on an open dense subset without being constant. 
This of course annihilates any hope of computing GCDs of polynomials when only approximations of 
them are known. Fortunately, there exists a standard way to recover continuity in this context: it 
consists in replacing GCDs by subresultants which are playing an analoguous role. For this reason, 
in what follows, we will exclusively consider the problem of computing subresultants. 

Definitions and notations 

We recall briefly basic definitions and results about resultants and subresultants. For a more 
complete exposition, we refer to [2, §4.2], [ 8 , §3.3] and [12, §4.1]. Let 21 be an arbitrary ring and 
let A and B be two polynomials with coefficients in 21. We pick in addition two integers d,A and ds 
greater than or equal to the degree of A and B respectively. We consider the Sylvester application: 

2l<d B [A] x 2l<d j4 [A] 2i <dA+dB [X} 

(U, V) i-a AU + BV 

1 We observe that these formulas can be rephrased as follows: the absolute (resp. relative) precision on the result 
of a sum or a substraction (resp. a product or a division) is the minimum of the absolute (resp. relative) precisions 
on . 
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where 2l< d [X] refers to the finite free 21-module of rank d consisting of polynomials over 21 of degree 
strictly less than d. The Sylvester matrix is the matrix of ip in the canonical ordered basis, which 
are 

0),..., (X, 0), (1, 0), (0, X dA ~ l ), ..., (0,1)) for the source 
and (X dA+dB_1 ,..., X, 1) for the target. 

The resultant of A and B (computed in degree dA,ds) is the determinant of the ip; we denote it 
by Res dA,dB (A, B). We observe that it vanishes if dA > deg A or ds > deg B. In what follows, we 
will freely drop the exponent dA,ds if dA and ds are the degrees of A and B respectively. Using 
Cramer formulae, we can build polynomials U dA ’ dB [A,B) £ 2l< dB [AT] and V dA,dB (A, B) £ 2l< dA [J 5 f] 
satisfying the two following conditions: 

i) their coefficients are, up to a sing, maximal minors of the Sylvester matrix, and 

ii) A-U dA ’ dB (A,B) + B-V dA ' dB (A,B) = Res dA ’ dB (A,B). 

These polynomials are called the cofactors of A and B (computed in degree dA, ds)- 

The subresultants are defined in the similar fashion. Given an integer j in the range [0, d) where 
d = nhn(d J 4 , ds), we consider the “truncated” Sylvester application: 

iP r . 2l< dB _ i [X]x2l <dA _ i [X] 2t <dA+dB _ i [X]/2t< j [X] 

(U, V ) AU + BV. 

Its determinant (in the canonical basis) is the j- th principal subresultant of A and B (computed 
in degree dA,ds )• Just as before, we can construct polynomials U dA ’ dB (A, B) £ 21 <( i B -j[X] and 
Vf A ’ dB (A,B) £ 2l< (U _ j [X] such that: 

i) their coefficients are, up to a sing, maximal minors of the Sylvester matrix 2 , and 

ii) A-U dA ’ dB {A 1 B) + B ■V J dA ’ dB (A,B) = detip j (mod 21 <? -[X]). 

We set R dA ' dB ( A , B) = A - U dA ’ dB ( A , B) + B ■ V dA,dB (A, B): it is the j-th subresultant of A and B 
(computed in degree c?a, ds)■ The above congruence implies that R dA ’ B (A, B) has degree at most 
j and that its coefficient of degree j is the j-th principal subresultant of A and B. As before, we 
freely drop the exponent dA,ds when dA and ds are equal to the degrees of A and B respectively. 
When j = 0, the application ipj is nothing but ip. Therefore, ResQ A,di3 (A, B) = Res dA,d,B (A, B) 
and, similarly, the cofactors agree: we have Uq A,<1b (A,B) = U dA,dB (A, B) and Vg A,dB (A, B) = 
V dA ’ dB (A,B). 

We recall the following very classical result. 

Theorem 2.4. We assume that 21 is a field. Let A and B be two polynomials with coefficients in 
21. Let j be the smallest integer such that Res j(A,B) does not vanish. Then ReSj(A, B) is a GCD 
of A and B. 

Since they are defined as determinants, subresultants behave well with respect to base change: if 
/ : 21 —?> 21' is a morphism of rings and A and B are polynomials over 21 then Res^ A ’ ds (/(A), f(B)) = 

f(Res dA ’ dB (A, B)) where f(A) and f(B) denotes the polynomials deduced from A and B respec¬ 
tively by applying / coefficient-wise. This property is sometimes referred to as the functoriality 
of subresultants. We emphasize that, when / is not injective, the relation Res j(f(A),f(B)) = 
f(ReSj(A, B )) does not hold in general since applying / may decrease the degree. Nevertheless, if 
dA and ds remained fixed, this issue cannot happen. 

2 Indeed, observe that the matrix of ipj is a submatrix of the Sylvester matrix. 
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The subresultant pseudo-remainder sequence 

When 21 is a domain, there exists a standard nice Euclide-like reinterpreation of subresultants, 
which provides in particular an efficient algorithm for computing them. Since it will play an 
important role in this paper, we take a few lines to recall it. 

This reinterpretation is based on the so-called subresultant pseudo-remainder sequence which 
is defined as follows. We pick A and B as above. Denoting by (P % Q ) the remainder in the 
Euclidean division of P by Q , we define two recursive sequences (Si) and (c,) as follows: 

A, So = B, c_i = 1 

(-■Sj) Ei+ 1 s l : L 1 1 c~ ei ■ (Si -1 % Si) for * > 0 (9) 

' c -" £i+1 for * > - 1 . 

Here rq = deg 5), £i = n, i+ i — m and s t is the leading coefficient of Si if i > 0 and s_i = 1 by 
convention. These sequences are finite and the above recurrence applies until S t has reached the 
value 0 . 

Proposition 2.5. With the above notations, we have: 

Res j(A,B) = Si ifj = iii- i-l 

= 0 if ni < j < ni-i - 1 

= (^r) 6 ’ 1 ■ s i ifj = rii 

for all i such that Si is defined. 

Remark 2.6. The Proposition 2.5 provides a formula for all subresultants. We note moreover 
that, in the common case where n*_i = rij — 1, the two formulas giving Res ni (A, B) agree. 

Mimicing ideas behind extended Euclide’s algorithm, one can define the “extended subresultant 
pseudo-remainder sequence” as well and obtains recursive formulae for cofactors at the same time. 

Important simplifications occur in the “normal” case, which is the case where all principal 
subresultants do not vanish. Under this additional assumption, one can prove that the degrees of 
the Si’s decrease by one at each step; in other words, deg Si = ds — i for all i. The sequence (Si) 
then stops at i = ds- Moreover, the efs and the cf s are now all “trivial”: we have e% = 1 and 
Ci = Si for all i. The recurrence formula then becomes: 

S I+1 = a? • s" 2 ! • (Si—i % Si) for i > 1 . 

and Proposition 2.5 now simply states that Rj = Sd B -j■ In other words, still assuming that all 
principal subresultants do not vanish, the sequence of subresultants obeys to the recurrence: 

R d +i = A, R d = B, Rj -1 = r 2 • rj+j ■ (Rj+i % Rj) (10) 

where rj is the leading coefficient of Rj for j < d and r d+ 1 = 1 by convention. Moreover, a similar 
recurrence exists for cofactors as well: 

u d+1 = 1, U d = 0, Uj- 1 =r 2 j-rfl 1 -(Uj + 1 -QjUj) (11) 

V d+ \ = 0, U d = 1 , Vj-x = r 2 j ■ rj£j ■ (V j+1 - QjVj) ( 12 ) 

where Qj is quotient in the Euclidean division of Rj+i by Rj. 

Proposition 2.5 of course yields an algorithm for computing subresultants. In the normal 
case and assuming further for simplicity that the input polynomials are monic of same degree, 
it is Algorithm 2, which uses the primitive prem for computing pseudo-remainders. We recall 
that the pseudo-remainder of the division of A by B is the polynomial prem(A, B) defined by 
prem(A, B) = lc (B) degB ~ degA+1 (A%B) where lc(R) denotes the leading coefficient of B. 

Unfortunately, while working over a complete discrete valuation field K , the stability of Algo¬ 
rithm 2 is as bad as that of standard Euclide algorithm. The use of Algorithm 2 is interesting 
because it avoids denominators (i.e. we always work over W instead K) but it does not improve 
the stability. 


S -1 = 

Si+i = 
Q+l — 
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Algorithm 2: Subresultant pseudo remainder sequence algorithm 
Input : Two polynomials A,Bg Kd[X] (given at finie precision) 
Output: The complete sequence of subresultants of A and B. 


Rd 4— B', Td <— 1 
Rd~i 4— B — A 

for j = (d— 1 ), (d — 2 ),..., 1 do 
r : j f— coefficient in Xi of Rj 
if rj = 0 then raise NotlmplementedError 


R 


■j -1 


prem(Rj + i,Rj)/r 


.2 
i +1 


7 return Rd- 1 ,..., Rq 


Example 2.7. Applying Algorithm 2 with the input ( A 1 B ) of Example 2.3, we obtain: 

R 4 = (29 + 0(2 5 ))A 4 + (14+0(2 5 ))A 3 + (5 + 0(2 5 ))A 2 + (l7 + 0(2 5 ))A+ (l7 + 0(2 5 )) 

R 3 = (4 + 0(2 5 )) A 3 + (13 + 0(2 5 )) A 2 + (4 + 0(2 5 )) A + (16 + 0(2 5 )) 

R 2 = (5 + 0(2 5 )) A 2 + (20 + 0(2 5 )) A + 0(2 5 ) 

Ri = (1 + 0(2))A + (1 + 0(2)) 

Ro = 1 + 0(2) 

We observe in particular that the absolute precision on R 0 is 1, although it should be at least 5 
since Rq is given by an integral polynomial expression in terms of the coefficients of A and B. 
We note moreover that the relative precision on Rq (which is 1 as well) is worse that the relative 
precision we got on Sy (which was 3) while executing Algorithm 1 (c/ Example 2.3). 


3 Unstability of Euclide-like algorithms 

In this section, we provide strong evidences for explaining the average loss of precision observed 
while executing Algorithm 2. Concretely, in §3.1 we establish 3 a lower bound on the losses of 
precision which depends on extra parameters, that are the valuations of the principal subresultants. 
The next subsections (§§3.2 and 3.3) aim at studying the behaviour of these valuations on random 
inputs; they thus have a strong probabilistic flavour. 

Remark 3.1. The locution Euclide-like algorithms (which appears in the title of the Section) 
refers to the family of algorithms computed GCDs or subresultants by means of successive Euclidean 
divisions. We believe that the stability of all algorithms in this family is comparable since we are 
precisely loosing precision while performing Euclidean divisions. Among all algorithms in this 
family, we chose to concentrale ourselves on Algorithm 2 because it is simpler due to the fact that 
it only manipulates polynomials with coefficients in W. Nevertheless, our method extends to many 
other Euclide-like algorithms including Algorithm 1 ; this extension is left as an exercice to the 
reader. 

3.1 A lower bound on losses of precision 

We consider two fixed polynomials A and B with coefficients in W whose coefficients are known 
with precision 0(n N ) for some positive integer N. For simplicity, we assume further that A and 
B are both monic and share the same degree d. For any integer j between 0 and d — 1, we denote 
by Rj the j-th subresultant of A and B. 

In this subsection, we estimate the loss of precision if we compute the Rj ’s using the recurrence 
(10). In what follows, we are going to use a flat precision model : this means that a polynomial 

3 in a model of precision which is slightly weaker that the usual one; we refer to §3.1 for a complete discussion 
about this. 







P(X) is internally represented as: 


P(X) = a,i.X l + 0(n N ) with a,i £ K and N £ Z. 

2=1 

Iii other words, we assume that the software we are using does not carry a precision data on each 
coefficient but only a unique precision data for the whole polynomial. Concretely this means that, 
after having computing a polynomial, the software truncates the precision on each coefficient to the 
smallest one. One can argue that this assumption is too strong (compared to usual implementations 
of p-adic numbers). Nevertheless, it defines a simplified framework where computations can be 
performed and experiments show that it rather well reflects the behaviour of the loss of precision 
in Euclide-like algorithms. 

Let Vj be the valuation of the principal j-th subresultant of A, B and Wj be the minimum of 
the valuations of the coefficients of Rj . We of course have Vj > Wj and we set Sj = Vj — Wj. 

Proposition 3.2. Let A and B as above. Either Algorithm 2 fails or it outputs the subresultants 
Rj’s at precision 0(ir Nj ) with: 


Nj < N + Vj .)_i — 2 • (Sj -|_i + Sjj-2 + • • • + Sd- 1). 


Proof. Using that Rj+i and Rj have the expected degrees, the remainder (Rj+i % Rj) is computed 
as follows: 

we set: S = Rj+i — i"j+i ■ rj 1 ■ Rj 

and we have: Rj+i % Rj = S — s ■ rj 1 ■ Rj 

where s is the coefficient of degree j of S. Let us first estimate the precision of S. Using (7)—(8), 
we find that the computed relation precision on • rj 1 ■ Rj is min(iV 7 _ t _i — Vj+i, Nj — Vj). The 
absolute precision of this value is then M = min(7Vj + i — Sj, Nj — Sj + Vj +1 — Vj). This quantity is 
also the precision of S since the other summand Rj+i is known with higher precision. Repeating 
the argument, we find that the precision of (Rj+i % Rj) is equal to mm(M — Sj, Nj — Sj +val(s) — Vj) 
and therefore is lower bounded by M — Sj < Nj — 28j + V 3 +\ — Vj. From this, we derive A r y -_i < 
Nj — 25j — Vj + \ + Vj and the proposition finally follows by summing up these inequalities. □ 

The difference N — N 0 = —Vi + 2 l a l° wer bound on the number of digits lost after 
having computed the resultant using the subresultant pseudo-remainder sequence algorithm. In 
the next subsection (cf Corollary 3.6), we shall see that V\ and all Sj 's are approximative^ equal 
to on average. The loss of precision then grows linearly with respect to d. This confirms the 
precision benchmarks shown in Figure 1. We emphasize one more time that this loss of precision 
is not intrinsic but an artefact of the algorithm we have used; indeed, one should not loose any 
precision when computing resultants because they are given by polynomial expressions. 

3.2 Behaviour on random inputs 

Proposition 3.2 gives an estimation of the loss of precision in Euclide-like algorithms in terms of 
the quantities Vj and Sj. It is nevertheless a priori not clear how large these numbers are. The 
aim of this paragraph is to compute their order of magnitude when A and B are picked randomly 
among the set of monic polynomials of degree d with coefficients in W. In what follows, we assume 
that the residue field k = W/nW is finite and we use the letter q to denote its cardinality. 

We endow W with its Haar measure. The set 12 of couples of monic polynomial of degree d with 
coefficients in W is canonically in bijection with W 2d and hence inherits the product measure. We 
consider Vj, Wj and Sj as random variables defined on 12. 

Theorem 3.3. We fix j £ {0,..., d — 1}. Let Xq, ■ ■ ■, X^-i be d pairwise independant discrete 
random variables with geometric law of parameter (1 — 9 _1 ), be. 

P [Xi = k] = (1 — g -1 ) • q~ k (with 0 < i < d and k £ N). 

Then Vj is distributed as the random variable 
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d 

Yj — ^ ' min(Xj_j, Xj_.j_|_i, • • •, -^F/’+i) 
i=0 

with Xi = +oo if i < 0 and Xj = 0 if i > d. 

Remark 3.4. The above Theorem does not say anything about the correlations between the 
Xj’ s. In particular, we emphasize that it is false that the tuple (Vd_i,..., Vo) is distributed 
as (Yd_i,.. .,Yo). For instance, one can prove that (Vd-\,Vd- 2 ) is distributed as (X, X' + 
min(X', [X/2])) where X and X' are two independant discrete random variables with geometric 
law of parameter (1 — g _1 ) and the notation [•] stands for the integer part function. In particular, 
we observe that {Vd-i, Vd- 1 ) ^ (2, 1 ) almost surely although the events {V,i-\ = 2} and {Vd -2 = 1} 
both occur with positive probability. 

Nonetheless, a consequence of Proposition 3.10 below is that the variables Vj = t{Vj=o} are 
mutually independant. 

Theorem 3.5. For all j £ {0,..., d — 1} and all m £ N, we have: 


1% > m] > (g q ^_ 1 V m 

The proof of these two theorems will be given in §3.3. We now derive some consequences. Let 
a denote the following permutation: 


2 •• 

d 

2 

- + 1 

2 ' 1 

1+2 

3 •• 

• d- 1 

d 

d -2 

2 • 

d+ 1 

2 

d+ 3 

2 

<2+5 

2 

3 • 

d 

d- 1 

CO 

1 


if 2 Id 


if 2 \d. 


In other words, a takes first the odd values in [1, d] in increasing order and then the even values in 
the same range in decreasing order. 

Corollary 3.6. For all j £ {0, ..., d — 1}, we have: 


d ~j 1 

(!) E [Vj] = E a(i) _ i > in Particular < E [Vj] < 

i —1 ^ 


g 

(q-iy 


( 2 ) 


q-i — 1 
q 0 + 1-1 


< E[5j] < E[V)] 


(3) a[Vj} 2 = E ^(0^1)2 '< in Particular ^ < a[Vj] < 

(4) P [Vj >m}< q -m+0(V^) 

(5) E[max(V q, ..., Vd-i)] < log ? d + 0(^/log q d) 


q+q +T 
(g-i ) 2 


Proof. By Theorem 3.3, we have E[Vj] = Y^i =with ^ = min(Xj_j,. .., Xj + f) (j is fixed 
during all the proof). Our conventions imply that Zi vanishes if z > d — j. On the contrary, if 
i < d—j , let us define r(l),..., r(d—j) as the numbers rr(l),..., a(d—j) sorted in increasing order. 
The random variable Zi is then the minimum of t(i) independant random variables with geometric 
distribution of parameter (1 — q~ 1 ) and thus its distribution is geometric of parameter (1 — 

Its expected value is then and the first formula follows. The inequality < E[V)] is clear 

because —- is the first summand in the expansion of E[V)]. The upper bound is derived as follows: 


< E 

2 — 0 


1 oo 

i=0 


qi _ qi -1 


g 

(g-1) 2 ' 
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The first inequality of claim (2) is obtained from the relation 


oo oo 

E [<y = E m ■ p & = = E p [^ ^ 

m =1 m= 1 

using the estimation of Theorem 3.5. The second inequality is clear because 5j < Vj. 
The variance of Vj is related to the covariance of Zf s thanks to the formula 

Var(Vj-) = Y, Cov(Zj,Zj/). 

<d—j 


Moreover, given X and X' two independant variables having geometric distribution of parameter 
(1 — a -1 ) and (1 — 6 _1 ) respectively, a direct computation gives: 


Cov(X, min(A', A'')) = x y 


Applying this to our setting, we get: 


Co v(Zi,Zi.) 


( q e(i,i') _ 1)2 


where e{i,i') = min(r(i),r(i / )) = r(min(i, *')). Summing up these contributions, we get the 
equality in (3). The inequalities are derived from this similarly to what we have done in (1). 

We now prove (4). Let (.Zi)i>o be a countable family of independant random variable having all 
geometric distribution of parameter (1 — q~ 1 ). We set Z = niin(Zi,..., Zi). Cleary Vj < Z 
and it is then enough to prove: 

P [Z > m] < 

We introduce the event E m formulated as follows: there exists a partition (mi, ... ,rni) of m such 
that Xj > mi for all i < t. Up to a measure-zero subset, E m contains the event {Z > m}. We 
obtain this way: 

i 

P [Z > m] < P [E m ] < En P[A'i > mj 

i=1 

where the latter sum runs over all partitions (mi,... ,me) of m. Replacing P[Xi > m*] by q~ mi , 
we get P[U m ] < p(m) ■ q~ m where p(m) denotes the number of partitions of m. By a famous 
formula [1], we know that logp(?n) is equivalent to 7 ry / 2 m/ 3 . In particular it is in q °(V™) an d ( 4 ) 
is proved. 

We now derive (5) by a standard argument. It follows from (4) that 


P[max(Vb,..., U d _i)] < d ■ q ~ m +^ 


for some constant c. Therefore: 

OO 

E[max(Vo,..., U d _i)] < Y min ^’ d ' q~ m+c ^)- 

m= 1 


Let mo denote the smallest index such that d ■ q m o+ c %A" 0; j. e . mo — c^/mo > log d. Solving 
the latest equation, we get mo = log ? +0(yj log g d). Moreover J2m=m 0 d q _m+c%/ ™ is bounded 
independantly of d. The result follows. □ 


3.3 Proof of Theorems 3.3 and 3.5 

During the proof, A and B will always refer to monic polynomials of degree d and Rj (resp. Uj 
and Vj) to their j-th subresultant (resp. their j-th cofactors). If P is a polynomial and n is a 
positive integer, we use the notation P[n] to refer to the coefficient of X n in P. We set rj = Rj[j}. 

Preliminaries on subresultants. We collect here various useful relations between subresultants 
and cofactors. During all these preliminaries, we work over an arbitrary base ring 21. 
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Proposition 3.7. The following relations hold: 


. Uj—\Vj - UjVj-i = (-1 Vr]; 

• Uj[d-j-l] = -Vj[d-j- 1] = (~l) j r j+1 ; 

• (Rj,Rj-i) = rf J ~ k ~ 1) R k for k < j; 

• Res = r ? d - j - k - 1) U d - t - k fork<d- j. 

Moreover rj depends only on the 2 (d — j) — 1 coefficients of highest degree of A and B. 

Proof. By functoriality of subresultants, we may assume that 21 = Z[ao,..., ad- i, bo,..., bd-i] and 
that A and B are the two generic monic polynomials A = X d +Y^i=o a-iX 1 and B = X rf +JT_ 0 biX 
Under this additional assumption, all principal subresultant are nonzero. Therefore, the sequences 
( Rj)j , ( Uj)j and {Vj)j are given by the recurrences (10)-(12). The two first announced relations 
follow easily. Let now focus on the third one. We set Rj = Rj and R k = r - R k for k < j. 

An easy decreasing induction on k shows that this sequence obeys to the recurrence: 

Rk -1 = ffc ■ ?k +1 • {Rk+l % Rk) 


where fj = 1 and f k is the coefhcient of R k of degree k for all k < j. Comparing with (10), 
this implies that R k is the fc-th subresultant of the pair (Rj, Rj—i) and we are done. The fourth 
equality and the last statement are proved in a similar fashion. □ 

For any fixed index j e {1,..., d— 1}, we consider the function ipj that takes a couple (A, B) e 
21 d[X } 2 to the quadruple (Uj, Uj-i, Rj, Rj-i). It follows from Proposition 3.7 that ipj takes its 
values in the subset £j of 

(21 < d -j-i[X]) x (21 <d-j[X]) x (21 <j(X}) x (21 <^[1]) 
consisting of the quadruples (U :j ,Uj-\, IZj, IZj _i) such that: 

Uj-i (d-j) = (-I?" 1 TIM 

and Res d ~ :i ’ d ~^~ 1 (Uj-i, Uj) = -Uj [j] 2 ^-!- 1 ). 

Let £* be the subset of £ 3 defined by requiring that 7 Zj [j] is invertible in 21. In the same way, we 
define fi* as the subset of 21 d[X ] 2 consisting of couples (A,B) whose j-th principal subresultants 
(in degree (d,d)) is invertible in 21. 

Proposition 3.8. The function ipj induces a bisection between f l* and £j. 

Proof. We are going to define the inverse of ipj. We fix a quadruple (U 3 , U 3 -\, 1Z 3 , Rj-i) in £* and 
set a = IZj [j]. Let W 3 and W ; _i denote the j-th cofactors of (it ; _i, U 3 ) in degree (d—j,d—j— 1). 
Define V 3 = aWj and Vj-i = — aW ; _i where a = a 4 l- 4d + 6 . The relation: 

Uj-iVj — UjVj-i = a 2 . (13) 

then holds. We now define A and B using the formulae: 

f A = (~iy ■ a -2 • (V^-r - Vj-jTlj) 

{ B = (-1 y - 1 ■ a ~ 2 ■ (UjUj-j - Uj- X lZj) 

and let ipj be the function mapping (Uj,Uj-i,TZj,lZj-i) to ( A,B ). The composite tpj oipj is easily 
checked to be the identity: indeed, if ipj(A,B) = (Uj,Uj-i,lZj,lZj-i), the relation (13) implies 
that Vj_i and V 3 are the missing cofactors and, consequently, A and B have to be given by the 
system (14). 

To conclude the proof, it remains to prove that the composite in the other direction ipj o ip 3 is 
the identity as well. Since both ipj and ipj are componant-wise given by polynomials, we can use 
functoriality and assume that 21 is the held Q(cq, c\,..., c n ) (with n = 2 d) and that each variable 
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Cj corresponds to one coefficient of U 3 , U 3 -\, IZj and TZj-i with the convention that Co (resp. 
(—l)' 7_1 co) is used for the leading coefficients of IZj (resp. U 3 -\). Set: 

(A, B) = ipj ( Uj , Uj ~\, IZj , Hj- 1 ) 
and (Uj, Uj-i, Rj, Rj-i) = ipj(A, B) 

Since 21 is a field and 7 Zj\j\ does not vanish, the Sylvester mapping 


2l< d _ j [A] x2l< d _,-[A] 2l< 2d _ j [A]/2l< j [A] 

( U , V) AU + BV 

has to be bijective. Therefore there must exist A G 21 such that IZj = X ■ Rj and Uj = A • U 3 . 
Similarly (JZj-\,Uj-{) = p ■ (Rj-i,Uj-±) for some /r G 21. Identifying the leadings coefficients, we 
get X = p. Writing Res d_J ’ d_J_1 (ZYj_i,ZY 7 ) = Res d_J ’ d_ - , ' _1 (?7 3 -_i, Uj), we get A 2 ^ - - 7 )-! = \ Since 
the exponent is odd, this implies A = 1 and we are done. □ 

Corollary 3.9. We assume that 21 = W. Then the map ipj '■ H* —► £* preserves the Haar 
measure. 

Proof. Proposition 3.8 applied with the quotient rings 21 = W/'K n W shows that (ipj mod 7r”) is a 
bijection for all n. This proves the Corollary. □ 

The distribution in the residue field. We assume in this paragraph that 21 is a finite field of 
cardinality q. We equip fig = 21 d[A] 2 with the uniform distribution. For j G {0,..., d — 1} and 
(A, B) G fist, we set Vj(A,B) = 1 if rj(A,B) vanishes and Vj(A,B ) = 0 otherwise. The functions 
Vj ’s define random variables over fig ■ 

Proposition 3.10. With the above notations, the Vj ’s are mutually independant and they all follow 
a Bernoulli distribution of parameter A. 

Proof. Given J C {0, ...,d — 1}, we denote by flg(J) the subset of fig consisting of couples 
(A, B) for which rj(A,B) does not vanish if and only if j G J. We want to prove that flg(J) 
has cardinality q 2d ~ CardJ (q _ i)Card j ^his. W e introduce several additional notations. 

First, we write J = {n\,... ,ng} with n\ > n 2 > • • • > ne. and set nt +1 = 0 by convention. 
Given n and m two integers with m < n, we let V( mjn ) denote the set of polynomials of the form 

a m X m +a m+ iA' m+1 -hanA'” with a* G 21 and a n ^ 0. Clearly, V( m) „) has cardinality (q— l)q n ~ m . 

If P is any polynomial of degree n and to < n is an integer, we further define P[m:\ G as 

the polynomial obtained from P by removing its monomials of degree < m. Finally, given (A,B) 
in fig, we denote by (Sj(A, B)) its subresultant pseudo-remainder sequence as defined in §2.2. We 
note that, if (A,B) G fig (J), the sequence (Si(A,B)) stops at i = i and we have deg Si = m for 
all i. We now claim that the mapping 

Aj '■ fla (J) V( ni,ri2) ^ ,ne+ 1 ) 

(A,B) (-)• (Si(A,B)[n.i+ 1 : ]) 

is injective. In order to establish the claim, we remark that the knowledge of Si-\(A,B ) and 
Si(A,B)[rii+ 1 :] (for some i) is enough to reconstruct the quotient of the Euclidean division of 
Si(A,B) by Si-\(A,B). Thus, one can reconstruct St(A,B) from the knowledge of Si- 2 (A,B), 
Si-i(A,B) and Si(A,B)[rii+ 1 :]. We deduce that A j(A,B) determines uniquely all Si(A,B)’ s and 
finally A and B themselves. This proves the claim. 

To conclude the proof, we note that the claim implies that the cardinality of flg(J) is at most 
q 2d ~ e (q — 1) £ . Summing up these inequalities over all possible J, we get Cardflg < q 2d . This latest 
inequality being an equality, we must have Card flg( J) = g 2d ~ Gard J (q — ijCard j f or a p j □ 

Proof of Theorem 3.5. We assume first that j < d— 1. Proposition 3.10 above ensures that r 3+ i 
is invertible in W with probability (1 — q -1 ). Moreover, assuming that this event holds, Corollary 
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3.9 implies that Rj is distributed in W<j[X] according to the Haar measure. An easy computation 
gives P[<5j > m | r J+ i £ W x ] = and therefore: 


P[<5j > to] > (1 


W 1 )- 


q(q J - 1) 

Y+ 1 -1 


(q- i)(g J -!) 
Y+ 1 - 1 


The case j = d — 1 is actually simpler. Indeed, the same argument works expect that we know for 
sure that i'j+i = rj, is invertible since it is equal to 1 by convention. In that case, the probability 
is then equal to ■ 

Proof of Theorem 3.3. We fix j £ {0,. .., d — 1}. We define the random variable Vj 0> as the 
greatest (nonnegative) integer v such that all principal subresultants ry have positive valuation for 
j' varying in the open range (j — v, j + v) (with the convention that ry = 0 whenever j' < 0). It 
is clear from the definition that rj— v or rj +v (with v = V^ 0 " 1 ) has valuation 0. Moreover, assuming 
first that val(r J+v ) = 0, we get by Proposition 3.7: 


val(rj) = v + val(rj ’ J+1 (A^°\l?( 0 ))) 
with = r . +v x3--»-i ■ R j+v\j-v-l :], 

and B (1) = ■ Rj +v -i[j-v-l :] 


where we recall that, given a polynomial P and an integer m, the notation P[m:\ refers to the 
polynomial obtained from P by removing its monomials of degree strictly less than m. We notice 
that all the coefficients of B^ lie in W because ry has positive valuation for j' £ (j — v,j + v). 
Furthermore, Corollary 3.9 shows that the couple (A 1 ), £?W) is distributed according to the Haar 
measure on (W 2 V -i[X]) 2 . If val(ryi_„) = 0, one can argue similarly by replacing Rj+ V and Rj + V -1 
by the cofactors Uj- V and Uj - V+ 1 respectively. Replacing (A , B) by (A 1 ), B^), we can now define 
a new random variable vj 1 ' J and, continuing this way, we construct an infinite sequence vj m ' > such 
that Vj = E m >o v j m) - 

We now introduce a double sequence (A’| m '*)o<i<d im >o of mutually independant random vari¬ 
ables with Bernoulli distribution of parameter ^ and we agree to set = 0 for j' < 0 and 

Xj™'* = 1 for j > d. It follows from Proposition 3.10 (applied with 21 = k) that has the 


same distribution than yJ 0 ' 1 = ^“ =1 min(XjY,..., Aj+'f). In the same way, keeping in mind 
that A^> and B^) have both degree 2 Vj°^ — 1, we find that Vj' 1 ' 1 has the same distribution than 


do) 


(o) ■ 


YaIi 1 min {X^% X^), which can be rewritten as y/ 1} = J2i=i min(AYY, AjY,..., X^X^). 

More precisely, the equidistribution of (A^ l \B^) shows that the joint distribution is 

the same as those of (Y^°\ Y)^). Repeating the argument, we see that is distributed 

as (Y (m) ) m >o where: 


do) xd 1 ) 


do) y (i)- 


Y j (m) = Y, min (Xj% 


v'( m ) y(0) 

^j-i ’ • • • > ■ j ^j+i )• 


Setting finally Xi = Yl m >o min(X}°\..., X^), we find the AYs (0 < i < d ) are mutually 
independant and that they all follow a geometric distribution of parameter (1 — q^ 1 )- We now 
conclude the proof by noting that Yj equals l min(X,_j,..., Xj +i ) (recall that the Xf m ^’s only 

take the values 0 and 1). 


4 A stabilized algorithm for computing subresultants 

We have seen in the previous sections that Euclide-like algorithm are unstable in practice. On the 
other hand, one can compute subresultants in a very stable way by evaluating the corresponding 
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Figure 3: Picture of a lattice in the ultrametric world 


minors of the Sylvester matrix. Doing so, we do not loose any significant digit. Of course, the 
downside is the rather bad efficiency. 

In this section, we design an algorithm which combines the two advantages: it has the same 
complexity than Euclide’s algorithm and it is very stable in the sense that it does not loose any 
significant digit. This algorithm is deduced from the subresultant pseudo-remainder sequence 
algorithm by applying a “stabilization process”, whose inspiration comes from [6] . 

4.1 Crash course on ultrametric precision 

In this subsection, we briefly report on and complete the results of [6] where the authors draw the 
lines of a general framework to handle a sharp (often optimal) track of ultrametric precision. In 
what follows, the letter W still refers to a complete DVR while the letter K is used for its fraction 
field. 

4.1.1 The notion of lattice 

As underlined in Remark 2.2, the usual way of tracking precision consists in replacing elements of 
W — which cannot fit entirely in the memory of a computer — by balls around them. Using this 
framework, a software manipulating d variables in W will work with d “independant” balls. The 
main proposal of [6] is to get rid of this “independance” and model precision using a unique object 
contained in a d-dimensional vector space. In order to be more precise, we need the following 
definition. 

Definition 4.1. A IV-lattice in a finite dimensional vector space E over AT is a lU-submodule of 
E generated by a A'-basis of E. 

Although the defintion of a lattice is similar to that of Z-lattice in K. d , the geometrical rep¬ 
resentation of it is quite different. Indeed, the elements of W themselves are not distributed as 
Z is in M but rather from a ball inside K (they are exactly elements of norm < 1). More gener¬ 
ally, assume that E is equipped with a ultrametric norm || • ||b compatible with that on K ( i.e. 
||Ax||e = |A| • ||x||_e for A € K, x £ E). (A typical example is E = K n equipped with the sup 
norm.) One checks that the balls 

B E {r) = { x G E | ||x||s < r } 

are all lattices in E. Moreover, any lattice is deduced from B e(1) by applying a bijective linear 
endomorphism of E. Therefore, lattices should be thought as special neighborhoods of 0 (see 
Figure 3). As a consequence, cosets of the form x + H, where H is a lattice, appear as interesting- 
candidates to model precision. This feeling is consolidated by the following result which roughly 
speaking claims that such cosets behave quite well under differentiable maps. 
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Lemma 4.2 ([6], Lemma 3.4). Let E and F be two normed finite dimensional K-vector spaces. 
Let f : E —»• F be a function of class C 1 and let x be a point in K n at which the differential of f, 
denoted by f'{x), is surjective. Then, for all p £ (0,1], there exists 5 > 0 such that the following 
equality holds: 

f(x + H) = f(x) + f(x)(H) (15) 

for any lattice H satisfying B E {pr) C H C B E (r) for some r < S. 

In what follows, we will often use Lemma 4.2 with p = 1. It states in this particular case that 

f[x + B E {r )) = f(x) + f'(x){B E {r )) (16) 

as soon as r is small enough. It is moreover possible to provide an explicit upper bound on r 
assuming that / has more regularity. The case of locally analytic functions is treated in [6] in full 
generality. Nevertheless, for the application we have in mind, it will be enough to restrict ourselves 
to the simpler case of integral polynomial functions. In order to proceed, we assume that E is 
endowed with distingushed “orthonormal” basis 4 , that is a basis (ei,...,e n ) with the property 
that II E"=l x i 6 i\\E = maxi<j<„|xi| for all families of A,;’s lying in K. In other words, the choice of 
this distingushed “orthonormal” basis defines a norm-preserving isomorphism between E and K n 
endowed with the sup norm. We assume similarly that we are given a distingushed “orthonormal” 
basis (/i,..., f m ) of F. Then any function / : E —> F can be written in our distinguished system 
of coordinates as follows: 


f(x) = ^2 F j{ x i> ■ ■ ■, x n )fj with x = ^2 x i e i- 
j =i *=i 

Definition 4.3. The function / is integral polynomial if all Fj’s are polynomials functions with 
coefficients in W. 

Example 4.4. Let us examine more closely the case of polynomial spaces since it will be considered 
repeadtly in the sequel. We take E = K <n [X] and F = K <m [X] and endow both with the Gauss 
norm, which is defined by: 


IIo-o + aiX H-1- a n -\X n 1 || E = max(|a 0 |, |ai|,..., |a„_i|) 

IIfro + hX 4-f b m - iX m ~ 1 || i r = max(|fr 0 |, |&i|,..., |fr m _i|) 

It is clear from these definitions that the canonical basis (1, A,... , X" -1 ) and (1,X,... ,X m_1 ) 
of E and F respectively are “orthonormal”. Moreover the coordinates in these basis are the af s 
and the bfs respectively. Hence, an integral polynomial function / : E —> F is nothing but a 
function mapping a a polynomial P to a polynomial Q whose coefficients are given by polynomial 
expressions which involve only the coefficients of P and some constants in W. 

Obviously, all integral polynomial functions are function of class C 1 (and even locally analytic), 
so that Lemma 4.2 applies to them. Proposition 4.5 below exhibits an exploit value for the bound 
5 appearing in Lemma 4.2 when / is integral polynomial and r = 1. 

Proposition 4.5. Let f : E —» F be an integral polynomial function and x £ B E (1). Then, 
Eq. (16) holds as soon as B E (r) C f(x)(B E (l)). 

Proof. It is a direct corollary of [6, Proposition 3.12]. □ 

4.1.2 Application to precision 

Let us now briefly explain how Lemma 4.2 can be utilized for tracking precision. 

4 One can prove that such a basis always exists. 
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Tracking precision locally Assume first that we want to perform a given rather simple op¬ 
eration - corresponding, say, to an elementary step ( e.g. an iteration of the main loop) of the 
algorithm we are executing — modeled by a function g of class C 1 defined on an open subset U 
of a finite dimensional normed A'-vector space E and taking values in another finite dimensional 
normed A'-vector space F. Our input is an approximated element of U which is represented by 
a coset C with respect to some lattice H, that is a subset of U of the form C = x + H for some 
x £ U. We would like to insist on the following: the value of a; is a priori not given; only is given 
the subset C. However, since H is stable under addtion, we have C = x + H for any clement 
x £ C. 5 As explained in §2.1.2, assuming that g is given as an algebraic expression, the naive 
solution for evaluating g[C) consists in using formulas (5)-(8). However, this often results in an 
overestimation on the precision, in the following sense: this method leads to some inclusion 

g(C) = g(x + H) C y + H naive 

where y £ F and H na i ve is a lattice which is generally much more larger that g'(x)(H), the latter- 
being the best possible one according to Lemma 4.2 (assuming that the assumptions of this Lemma 
are fullfiled). In order to avoid this and be sharp on precision, another solution consists in splitting 
the computation of g(C ) into two parts as follows: 

(A) compute g'(x)(H), and 

(B) compute g(x) for some x £ C. 

Part (A) is not easy to handle in full generality: in order to be efficient, a special close analysis 
taking advantage of the particular problem under consideration is often necessary. For now, let us 
simply assume that we have given two lattices H m j n and A max with the property that: 

H min Cg\x)(H)cH max . (17) 

We shall see later (c/ §4.2) how these lattices can be constructed — for a negligible cost — in the 
special case of subresultants. 

We now focus on part (B), which also requires some discussion. Indeed, computing g(x) is 
not straightforward because x itself lies in a A'-vector space and therefore cannot be stored and 
manipulated on a computer. Nevertheless, one can take advantage of the fact that x may be chosen 
arbitrarily in C. More precisely, we pick a sublattice H' of H and consider the new approximated 
element x + H' C x + H. Concretely, this means that we arbitrarily increase the precision on the 
given input x. Now, applying the naive method with x + H ', we compute some y £ F and some 
lattice H ' naive C F with the property that: 

g(x + H') C y + H ' naive . 

If furthemore H' is chosen in such a way that H aaive C H min , the two cosets y + g'(x)(H) and 
g{C) have a non-empty intersection because g{x) lies in both. Therefore they must coincide. We 
deduce that y £ g{C). This exactly means that y is an acceptable value for g(x) and we are done. 
Moreover, estimating the dependance of H' naive in terms of H' is usually rather easy (remember 
that g is supposed to model a simple operation). Hence since H mm is known — as we had assumed 
— finding H' satisfying the required assumption is generally not difficult. 

Tracking precision globally As already said, we shall use the above method for tracking pre¬ 
cision while executing a single step in a complete algorithm. Let us now address the problem of 
“glueing”. We consider an algorithm F consisting in a succession of n steps G 0 ,..., G n -i- It is 
modeled by a function / : U —> F of class C 1 where U is an open subset in a finite dimensional 
normed A-vector space E and A is a finite dimensional normed A'-vector space. The input of F is 
an approximated element in U represented as a coset C = x + H where x £ U and H is a lattice. 
We also introduce notations for each individual step. For all i, we assume that G; is modeled by a 

5 This assertion means that any element of the “rectangle” C is a center of it... which might be surprising if we 
are accustomed to real numbers. 
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Figure 4: Method for tracking precision based on Lemma 4.2 


function gi : Ui -» {/*+i of class C 1 where Ui is an open subset is some normed A'-vector space A, 
and, by convention, Uq = U, E 0 = E and U n = E n = F. We thus have: 

/ = g n - 1 o g n _ 2 o---og 1 og 0 . 

For all i, we set fi = gi_\ o • • • o g 0 . It is the function modeling the execution of the i first steps 
of our algorithm. We further define Xi = fi{x) and Hi = f[(x){H). The chain rule for composing 
differentials readily implies the recurrence 

H i+ 1 = S /' i (x i )(H i ) (18) 

For simplicity, we make the following assumptions: 

• the Zp-submodule Hi is a lattice in Ei such that Xi + Hi C [/*; 

• the triple (^, Xi, Hi) satisfies the assumptions of Lemma 4.2; 

• for all i, we have succeeded in finding (good enough) explicit lattices A m i n ,i and H maXt i such 
that H m - ln i C Hi C H max i ; 

• for all i, we have succeeded in finding an explicit lattice H[ such that, while tracking naively 
precision, we end up with an inclusion 

4“ H^) Xi-\- 1 H na { ve i-i-\ 

with H na [ ve i^i C H m i n ,z+i- 

We note that the first and the second assumptions are quite strong because they imply in particular 
that the sequence of dim Ei is non-increasing. However, it really simplifies the forthcoming discus¬ 
sion and will be harmless for the application developed in this paper. As already mentionned, the 
construction of and A maXi , will generally follow from a theoretical argument depending on 

the setting, while exhibiting H[ will often be straightforward. Anyway, we are now in position to 
apply the method for tracking precision locally we have discussed earlier to all g^s. This leads to 
a stabilized, version of the algorithm F whose skeleton is depicted in Algorithm 3. 

The correctness of Algorithm 3 (under the assumptions listed above) is clear after Lemma 15. 

4.2 Application to subresultants 

We now apply the theory presented in §4.1 above to the problem of computing subresultants, i.e. 
the abstract Algorithm F is now instantiated to Algorithm 2. We split this algorithm into steps in 
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Algorithm 3: Stabilized version of F 
Input : x given at precision O(H) 
Output: g(x) given at precision 0(H maX}n ) 


1 

2 

3 

4 

5 


Xo <— X 

for i = 0 , ..., n — 1 do 

lift Xi to precision 0{H[) 
X%-\- 1 t G i(Xi) 
return x n + 0(Hmax,n) 


the obvious manner, each step corresponding to an iteration of the main loop. We thus consider 
the functions: 


g d : K d [X\ x K d [X\ -»■ 
(A, B) ^ 


K d [X] x A< d _ x [A] 
(■ B,A-B) 


and g-j : K< j+ 1 [X\ x K<j[X] K<j[X] x K<j- i[A] 

(Rj+i,Rj) t —> (Rj,Rj- 1) 

where f?j-i is defined as usual by Rj~\ = • r“ ) _ 1 • (-Rj+i %Rj) where rj (resp. Tj+fl) stands 

for the coefficient of degree j in Rj (resp. of degree j + 1 in Rj+\). We remark that gj is only 
defined on the subset consisting of pairs (Rj+i , Rj) for which Rj+i has degree j + 1; this reflects 
the fact that Algorithm 2 fails on inputs for which at least one principal subresultant vanishes. 
The composite function / = gi o ■ ■ ■ o g d (be careful with the order of the indices) models (a slight 
variant of) Algorithm 2. For all j, we put fj = gj+\ o • • • o g d ; it is the function: 

fj : K d [X] x K d [X] -»■ K<j[X\ x K^^X] 

(A,B) i->- (Res :) (A,f3),ReSj_i(A, B)). 

For simplicity, we assume in addition that the precision on the input ( A , B) is flat, meaning 
that all coefficients of A and B are known with the same absolute precision N. In the language 
of §4.1, this flat precision corresponds to the lattice H = tt n C where £ = W <( j[A] x W<d[A] 
is the unit ball in I\ d [X] x I\ d [X] with respect to the Gauss norm (c/ Example 4.4). Following 
§4.1, our first task consists in finding two lattices H mm j and H max j having the property that 

UnnnJ C /'(/t, «)(//) C // lmlx .,. 

Lemma 4.6. For all ( A,B ) £ K d [X] 2 , we have: 

r j ' A' C fj{A, B){C) c Cj 

where rj is the j-th principal subresultant of (A, B ) and Cj = W<j [A] x W<j- 1 [X] is the unit ball 
in K<j[X] x K<j_ i[A], 

Proof. The second inclusion is clear because fj is a polynomial function. Let us prove the first 
inclusion. One may of course assume that rj does not vanish, otherwise there is nothing to prove. 
Now, we remark that fj factors through the function 'ipj introduced in §3.3. By continuity, the j-th 
principal subresultant function does not vanish on a neighborhood of ( A,B ). By Proposition 3.8, 
tpj is injective on this neighborhood. Therefore so is fj. Furthermore, a close look at the proof of 
Proposition 3.8 indicates that a left inverse of fj is the function mapping (Sj,Sj- 1 ) to 

(-1 y ■ rf ■ (VjSj-t-Vj^Sj, -UjSj i +Uj-iSj) 

where Uj, Vj (resp. Uj- 1 , V)_i) are the j-th (resp (j — l)-th) cofactors of (A, B). Differenting this, 
we get the announced result. □ 

Lemma 4.6 ensures that one can safely take H mm j = r| • n N Cj and H max j = n N Cj. It finally 
remains to construct the lattice HI C K<j[X] x K<j- i[A]. For this, we remark that a naive 
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track of precision leads to a loss of at most 2 • val(r J+ i) digits while executing the step Gj (see also 
proof of Proposition 3.2 for similar considerations). Therefore, one can take H) = r 2 r 2 +1 ■ 7 t n C j. 
Instantiating Algorithm 3 in this particular case, we end up with Algorithm 4 below which then 
appears as a stable version of Algorithm 2. 


Algorithm 4: Stabilized version of Algorithm 2 
Input : Two polynomials A,Bg Kd[X] given at flat precision O^tt 11 ) 

Output: The sequence of subresultants of A and B given at flat precision 0{ 77 ”) 

1 Rd B', Td <— 1 

2 Rd— 1 B — A 

3 for j = (d — 1), (d — 2),..., 1 do 

4 r 3 -t— coefficient in X J of Rj 

5 if Vj > ^7 then raise NotlmplementedError 

6 lift (R j+1 ,Rj) at precision 0(n N + 2 v a i(r 3 -)+ 2 vai(r j+1 ))) 

7 |_ Rj^i prem(i?j + i, Rj)/Vj +1 

8 return Rd -1 + 0 ( 77 ^),..., Rq + 0(77^) 


Proposition 4.7. Algorithm 4 computes all subresultants of (A, B) at precision 0 {tt n ) under the 
following assumption 6 

(H); all principal subresultants of (A,B) do not vanish modulo tt n / 2 . 

R runs in 0(d 2 ■ M(A + max(Vo, • ■ ■, Vd- 1 )) bit operations where Vj denotes the valuation of r 3 and 
M(n) is the number of bit operations needed to perform an arithmetic operation (addition, product, 
division) in W at precision 0{ 77 "). 

Remark 4.8. In all usual examples (p-adic numbers, Laurent series), one can choose M(n) to be 
quasi-linear in n and the size of the residue field k. 

Proof. Correctness has been already proved (the assumption (H) ensures that Proposition 4.5 
applies to each <?.,■). As usual Euclide’s algorithm, Algorithm 1 requires 0(d 2 ) operations in the 
base ring W. Moreover, we observe that the maximal precision at which we are computing is upper 
bounded by N + 2max(Vo,..., Vd- 1 ). This justifies the announced complexity. □ 

According to Corollary 3.6, the expected value of the variable max(Vo,..., Vd-i) is in 0( log p d). 
Thus, the average complexity of Algorithm 1 is 0(d 2 ■ M(TV + logd)) bit operations. In all usual 
cases (cf Remark 4.8), this complexity is also 0(d 2 N ■ log \k\) bit operations. 

To conclude with, let us comment on briefly the hypothesis (H). We first remark that it is 
satisfied with high probability if N is large compared to 2 • log d p. Thus, replacing eventually N by 
3 • log d p (which does not affect the complexity), the assumption (b) is harmless on average - but 
maybe not on particularly bad instances. We moreover underline that, if we are just interested in 
computing the j- th subresultant for a particular j, then we just need to assume the non-vanishing 
of the principal subresultants in the range [j + 1, d — 1]. 

Open questions 

The first hypothesis we would like to relax is of course (H). Actually, it seems quite plausible that 
one can produce a stabilized version of the “complete” 7 subresultant pseudo-remainder sequence 
algorithm following the same strategy. Nevertheless, this extension is not completely straightfor¬ 
ward because designing it requires to understand precisely how the coefficients cf s (appearing in 
Eq. 9) alter the behaviour of the precision. We therefore let it as an open question. 

As it was presented, Algorithm 4 only accepts inputs consisting of a pair of monic polynomials 
having the same degree. It is actually not difficult to make it work with all couples of polynomials 

6 If this assumption is not fullfiled, the algorithms fails and returns an error. 

7 1.e. dealing with abnormal sequences as well. 
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(. A , B) such that lc(£?) is invertible in W and deg A > deg B. Indeed, it is enough for this to replace 
line 2 by: 

Rd- 1 «- (-l)degA-deg B (A%B). 

However, writing an extension of Algorithm 4 that accepts all inputs seems much more tricky and 
this is the second open question we would like to point out. 

Beyond this, one may wonder if one can use similar technics to compute not only subresultants 
but cofactors as well. For those indexes j such that r ; is invertible in W, the same analysis applies 
almost verbatim. However for other indexes j, the differential computation seems to be much more 
subtle. One can get around this issue by using lifting technics only when Tj is a unit in W and 
tracking precision naively otherwise: it is possible to get this way a stable algorithm whose average 
running time is acceptable but which seems to be bad in the worst case. Can we do better? 

Another quite interesting question is those of designing an algorithm which combines the preci¬ 
sion technology developed in this paper with the “half-GCD” methods. It is actually closely related 
to the previous question because “half-GCD” methods make an intensive use of cofactors in order 
to speed up the computation. 


5 Conclusion: towards p-adic floats 

When computing with real numbers, computers very often use floating point arithmetic. The rough 
idea of this model consists in representating all real numbers using the same number of digits (the 
so-called precision) and to apply rounding heuristics when final digits are unsettled. In comparison 
with arithmetic interval, floating point arithmetic has two main advantages. First, it allows simple 
and fast implementations. Second, experiments show that the obtained results have generally more 
much correct digits that those predicted by arithmetic interval. The counterpart is that, expect 
on small examples, obtaining proved results is generally intractable. 

In the p-adic setting, the analogue of floating point arithmetic has not been developed yet. One 
reason for this is probably the well-known saying: “in the p-adic world, rounding errors do not 
accumulate”. Consequently one might expect that interval arithmetic would provide sharp results. 
Nonetheless this hope is failing and examples are basic and numerous: p-adic differential equations 
[4, 10], LU factorization [7], SOMOS 4 sequence [6], resultants (this paper), etc. Consequently, 
interval arithmetic is not as good as one might have expected at first. Therefore, it probably makes 
sense to seriously study the analogue of floating point arithmetic in a ultrametric context. 

Let us describe quickly what might be this analogue and what are its advantages and disad¬ 
vantages. We keep the notations of the previous sections: the letter W denotes a complete discrete 
valuation ring with uniformizer tt and K is its fraction field. In the model of ultrametric float¬ 
ing point arithmetic, we fix a positive integer N (the precision ) and represent elements of K by 
approximations of the form: 

JV-l 

7 T e ■ XiW 1 (19) 

2=0 

where e is a relative integer and the Xi’s are elements of a fixed set of representatives of W modulo 
n with the convention that the representative of 0 € k is 0 £ W. We further assume that Xq / 0, 
i.e. e is the valuation of the sum (19). We see that this framework is quite similar to usual floating 
point arithmetics: the integer e plays the role of exponent, the uniformizer 7r plays the role of the 
basis and the value Y^iLo 1 plays the role of the significand (the mantissa). It remains to define 
operations © and 0 on approximations modeling addition and multiplication on K respectively. 
We do this as follows: given x and y two elements of K of the form Eq. (19), we compute x + y 
(resp. xy) in I\ , expand it as a convergent series Y2rL v Si7r * (with s v ^ 0) and define x © y (resp. 
x © y) by truncating the series at i = v + N. 

Similarly to real floating point arithmetic, the main advantages of ultrametric floating point 
arithmetic are the simplicity and the efficiency while the counterpart is the difficulty to get proved 
results. Moreover, the aforementioned examples are evidences that ultrametric floating point arith¬ 
metic may often compute much more correct digits than those predicted by an analysis based on 
interval arithmetic. In order to illustrate this last assertion, let us go back to the case of resultants 
discussed earlier in this paper. Let A and B be two monic polynomials of degree d (picked at 
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random) whose coefficients are all known at precision 0(n N ). We have proved that if we are using 
the model of interval arithmetic, then the subresultant pseudo-remainder sequence algorithm will 
output Res (.A, B ) at precision 0(ir N ~ Nint ) where N- Int grows linearly with respect to d on average. 
On the other hand, if we are using ultrametric floating point arithmetic, then the same algorithm 
will output Res (A,B) at precision 0(■jr N ~ Naoa - t ) where Vfl oat grows linearly with respect to logd 
on average. We emphasize furthermore that this result is proved\ From this point of view, floating 
point arithmetics seems to behave better in the ultrametric setting: we may hope to get proved 
results relatively cheaply. 
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